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described.  The  present  model,  a  modified  version  of  the  original 
model  by  Barnett  et  al .  (contract  N62306-68-C-0285,  U.S.  Naval 
Oceanog.  Off.,  1969)  includes  an  improved  representation  for 
weak,  nonlinear,  wave-wave  interactions  that  depends  on  spectral 
shape  and  prescribes  atmospheric  input  that  is  consistent  with 
measurements  by  Snyder  et  al .  (J.  FI.  Mech.,  102,  1-59,  1981). 

The  limiting  form  of  the  spectrum  is  defined  as  a 
Pierson-Moskowi tz  spectrum  with  a  variable  saturation  range 
parameter.  Some  additional  constraints  are  also  used  to  obtain 
stable  spectra  because  there  is,  in  effect,  a  mismatch  between 
the  numerous  degrees  of  freedom  of  the  two-dimensional  wave 
spectrum  and  the  level  of  sophistication  used  to  represent  the 
physical  processes  affecting  wave  growth. 

The  focus  of  the  present  note  is  on  the  wave  model  computer 
program.  The  note  is  organized  for  someone  wanting  to  implement 
the  program.  Hence,  top-level  flow  charts,  user  instructions, 
and  tape  output  formats  are  given.  In  addition,  a  detailed 
description  of  the  numerical  method  is  included. 
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A  DISCRETE  NONLINEAR  SPECTRAL  MODEL  FOR  OCEAN  WAVE  PREDICTION, 
DESCRIPTION  OF  COMPUTER  PROGRAM 


I.  INTRODUCTION 

Predicting  the  ocean  surface  wave  spectrum  is  a  formidable  task.  First,  a 
knowledge  of  the  temporally-  and  spatially -varying  winds  over  the  domain  of  interest 
is  required.  Second,  a  model  that  uses  these  winds  to  predict  the  wave  spectrum, 
which  depends  on  space,  time,  frequency,  and  direction,  is  needed.  Finally, 
computer  graphics  to  display  model  results  are  highly  desirable. 

For  the  present  note  the  winds  are  regarded  as  known.  Our  primary  purpose  is  to 
describe  in  considerable  detail  the  computer  program  that  was  developed  in  the 
present  study.  The  original  model  was  described  by  Barnett  (1968)  and  Barnett  et 
al.  (1969).  The  rationale  for  the  overall  model  is  given  in  Allender  et  al .  (1982). 
A  larger  perspective  on  wave  prediction,  as  well  as  intercomparisons  among  the 
present  and  other  models  in  ideal  situations,  can  be  found  in  Hasselmann  et  al . 
(1982).  The  remainder  of  this  note  is  organized  for  someone  who  wants  to  use  the 
model  program,  not  someone  who  wants  to  learn  about  wave  prediction  in  general.  The 
unfamiliar  reader  should  therefore  consult  the  preceding  references  for  the 
necessary  background  information.  In  the  final  section  of  the  note  we  take  a 
somewhat  broader  perspective  by  describing  some  of  the  pitfalls  we  found,  future 
changes  for  this  model,  and  possible  changes  for  wave  prediction  models  in  general. 

II.  DESCRIPTION  OF  COMPUTER  PROGRAM 
A.  WAVESET 

WAVESET  defines  the  problem  geometry  -  the  time  step,  spatial  grid, 
discrete  frequencies  and  directions,  boundaries,  and  ray  families.  These  quantities 


Influence  the  rays  In  two  ways:  (1)  In  principle  it  Is  possible  to  have  a  different 
family  of  parallel  rays  at  each  direction  for  each  frequency  (different  in  terms  of 
spacing  and  arrangement),  and  (2)  spacing  of  points  along  the  rays  is  based  on  the 
distance  CgAt,  where  Cg  is  the  group  velocity,  a  function  of  frequency  ,  and  *t 
is  the  time  step.  Ray  point  spacing  is  an  integral  multiple  (possibly  1)  of  this 
distance.  The  multiplier  can  be  specified  independently  for  each  frequency.  A 
different  set  of  ray  directions  can  be  specified  for  each  frequency. 

In  accordance  with  the  characteristic  equations  of  the  basic  partial 
differential  equation,  "(f,Q)  particles"  move  along  the  rays  (characteristics) 
from  ray  point  to  ray  point,  and  grow  in  a  certain  way.  The  arrays  of  these 
spectral  words,  and  other  arrays  (e.g.,  various  moments  of  the  spectra),  are 
initialized  in  WAVESET  and  stored  in  "tape"  (usually  disk)  files.  Tables  1A-B  show 
block  diagrams  of  WAVESET. 

B.  PROPAGR 

The  basic  function  of  PROPAGR  is  to  modify  the  spectral  density  of  and  on 
appropriate  time  steps  to  propagate  the  spectral  words.  (The  "energy"  carried  by 
these  particles  (words)  has  units  of  (length)2  (time)/radian. )  PROPAGR  consists 
mostly  of  a  set  of  nested  loops.  The  outermost  loop  is  a  time  step  loop.  One  pass 
through  this  loop  calculates  everything  that  takes  place  in  the  simulated  wave  field 
during  one  time  step  or  cycle. 

Within  the  time  step  loop  are  calculations  of  certain  quantities  that  are 
constant  for  that  time  step,  and  a  loop  over  all  frequencies.  At  the  end  of  each 
time  step,  when  the  directional  spectrum  has  been  evaluated,  it  is  possible  to 
calculate  spectral  moments,  such  as  fm*  and  Emax. 

♦Frequency  of  spectral  peak;  arbitrarily  initialized  to  0.25  Hz. 
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TABLE  1A 


FUNCTIONAL  BLOCK  DIAGRAM,  WAVESET 


Begin 

] 

Copy  input  cards  to  unit  2. 

General  initialization. 

Define  output^unit  LU  =  1.  $SET 

JO  =  0  •*—  99 

\$BOND 

Define  boundary  and  grid.  |$GRID 

Construct  geometry  record  & 
tables  for  neareslj  neighbor  search. 

Construct  grid  record  &  write  on  LU. 

Construct  boundary  record  l  write  on  LU. 

Restore  boundary  arrays. 

Write  geometry  record  on  LU. 

I 

Restore  grid  arrays. 
Frequency-direction  loop.* 

Construct  summary  record  &  write  on  LU. 

End  fi  e  LU. 

_ Were  there  _ _ 


I/O  errors 
? 


Stop  with 
message. 


♦Next  page. 
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TABLE  IB 


FREQUENCY-DIRECTION  LOOP,  WAVESET 


Loop  entrance 


*J  =  frequency  index,  K  =  direction  index. 
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TABLE  2A 


FUNCTIONAL  BLOCK  DIAGRAM,  PROPAGR 


C  Begin  ) 


Initialize. 

Copy  input  cards 

to  unit  2. 

Copy  WAVESET  output  tape 

to  unit  11. 

Construct  initial 

summary  tape 

(unit  21). 

Define: 

LIN  =  11 

LOUT  =  12 

JIN  =  21 

JOUT  =  22 

_ * _ 

Time  step  loop.* 


C  Stop*  ) 


$DBTA 


*Next  page. 

+Stop  test  &  stop  contained  in  SUBROUTINE  ADVANCE  (in  time  step 
loop) . 


TABLE  2B 

TIME  STEP  LOOP,  PROPAGR 


Loop  entrance 


Set  stop  indicator  off.* 


/ Stop  X  y 
indicator^, 
v  on?  / 


Loop  exit 


'Copy  grid  record  and  boundary  record  to  LOUT. 


Wind  input. 


Copy  geometry  record  to  LOUT. 


Process  summary  record. 

Calculate  y,c<  ,  fw  for  each  station. 
Convert  W  to  degrees. 


Initialize  spectrum  output  routine  for  this  step.  I  $SPRT 


Frequency  loop.+ 


Generate  new  summary  record  &  write  on  LOUT  &  JOUT. 
Calculate  moments,  fm,  Emax  for  each  station. 


Rewind  files. 

Exchange  identities:  LIN^LOUT,  JIN^JOUT. 


SUBROUTINE  ADVANCE* 

Increment  time. 

Produce  summary  printout. 

Turn  stop  indicator  on  if  stop  time  attained. 


*Stop  test  &  stop  contained  in  SUBROUTINE  ADVANCE 
+Next  page. 


TABLE  2C 


FREQUENCY  LOOP,  PROPAGR 


♦Next  page. 
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TABLE  2D 


DIRECTION  LOOP,  PROPAGR 
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Within  the  frequency  loop  there  are  calculations  of  quantities  constant  for  a 
given  frequency  and  there  is  a  loop  over  directions.  The  calculations  that  directly 
Involve  the  spectral  words  take  place  within  the  direction  loop.  Propagation  and 
interaction  take  place,  and  the  summations  for  the  Id  spectram  and  the  moments  are 
accumulated. 

Two  types  of  files  (external  storage)  are  used  to  pass  the  complete  state  of 
the  system  from  one  cycle  to  the  next.  The  cycle  tape  contains  information  on  the 
grid,  boundaries,  details  of  the  geometry,  quantities  constant  for  each  frequency, 
spectral  words  and  related  information,  integrals  of  energy  over  direction,  and 
integrals  over  both  direction  and  frequency.  The  spectral  words  and  integrals  are 
also  written  on  a  summary  tape.  One  pair  of  these  files  serves  as  input  to  a  cycle; 
another  pair  is  output.  In  addition,  for  output  purposes,  a  spectral  summary  tape 
is  created.  This  file  contains  the  directional  spectrum  at  each  grid  point  (found 
by  nearest  neighbor  interpolation).  The  subroutine  that  writes  this  file  also  can 
print  the  Id  and  2d  spectra  at  selected  grid  points.  Tables  2A-D  show  block 
diagrams  of  PROPAGR. 

III.  CALCULATION  OF  SOURCE  FUNCTIONS  FOR  SPECTRAL  DENSITY 
A.  MOMENTS 

The  following  quantities  are  calculated  at  each  grid  point  (using  nearest 
neighbor  interpolation  on  ray  point  values  where  needed). 

1.  The  variance  (energy)  EAV,  which  includes  an  f-5  tail,  is  defined  by 

Esum.  -  21  c!L  F  (  £,  9)  2-tt  a 

•j  k 

.OOfrl 

i,c  vun*  l 

EAV  =  ESUM  +  SE 
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Here  fmax  Is  the  highest  discrete  frequency  used  in  the  model,  not  the  frequency 
of  the  spectral  peak,  and  g  is  the  acceleration  of  gravity. 

2.  The  mean  frequency  FQAV  accounting  for  an  f-5  tail  is  found  from 

FglSiam  -  5  ^  £  &-T  Xlt  £>  6^ 

j  fc 

q  z  *5"  _ 

C"-t  t2'"' )  L"^>AV  ^  C  2  ~  ^  /lOO  J  ^ 

FQAV  =  (FQSUM  +  SF0)/EAV 

3.  The  mean  wave  direction  THAV  is 

ThSmMS  =  2T  ZL  P  (  ^ 

-j  K 

TliiiAMC  r  PH, 8)  tos  kfj 

J  fa. 

THAV  =  tan-1  (THSUMS/THSUMC) 

4.  The  rms  directional  spread  S  is 

Ti+R>u./^-  2,T6dA 

J  fa. 

S  =  THRMS  =  (THRSUM/EAV)1/2 
and  0  =  THAV. 

B.  OTHER  SPECTRAL  PARAMETERS 

Additional  parameters  needed  to  evaluate  the  source  functions  are  found 


2ir  &dk 


2-Tr  4 


as  follows. 


1.  The  peak  frequency  of  the  1-D  spectrum  FM  and  the  maximum  E^x  are 
found  by  fitting  a  quadratic  to  the  discrete  peak  and  evaluating  the  quadratic  at 
Its  maximum.  (Assuming  the  discrete  peak  occurs  at  subscript  2  (arbitrary  choice) 
the  formulas  are: 

*0  =  CF i(f if32  -  f 3^22 )  +  ^2^ ^3f i2  -  fif32) 

+  F3 { f l f 22  -  f2Fi2)3/D 
ai  =  CF 2f 32  -  F3f22)  +  (F3F12  -  FXf32) 

+  ( F 1 f 22  -  F2fi2)]/D 
a2  =  C ( f 2^ 3  -  f 3F 2 )  +  (f 3F1  -  f 1F3 ) 


+  ( f iF2  -  f 2F ! )]/D 

1  f!  f^ 

D  =  det  1  f2  f2* 

1  f3  f3J 

The  derived  location  and  value  of  the  maximum  then  are 
Fm  =  -ai/(2a2) 

EMAX  =  a0  -  ai2/(4a2) 

2.  The  saturation  range  (or  energy  level)  parameter  <<  for  use  in  the  limiting 
form  of  the  spectrum,  the  nondimensional  peak  frequency  and  the  wind  frequency 
fu  are  defined  as: 

V  =  ANU  =  U  FM/g 


<*  =  ALP  =  0.032  V2/3,  -0081  _<  •<  5  .02 
^  (FW  «  0.13  g/(/,  O  *  0 

0  (fw  =  1.0  ,  u  =  0 
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3.  The  peakedness  GAMA  of  the  1-d  spectrum  compared  to  a  Pierson-Muskowitz 
spectrum  EPM  is 

EPM  =  176.5  ALP/(FM)5 
=  GAMA  =  EMAX/EPM,  1.0  <  GAMA  <  4.0 

4.  The  scaling  term  to  convert  refrence  values  of  the  nonlinear  transfer 
Sni  to  values  appropriate  for  an  evolving  model  spectrum  is 

A*-**- 

C  »  f-  0.7.3S-  (  GftMA  -  /  )  ] 

C.  ATMOSPHERIC  INPUT 

The  form  of  S^n  used  in  the  oresent  model  yields  an  exponential  growth 
mechanism.  The  linear  growth  mechanism,  operative  in  the  old  model,  is  no  longer 
used.  The  source  function  has  the  same  form  as  described  by  Barnett  (1968),  but  has 
been  reduced  by  a  factor  of  four;  the  nonlinear  transfer  (Sec  III -D)  makes  up  the 
difference  in  the  present  model.  The  source  function  is  simply: 

Sin  =  0.25  (iY 

where  Y  is  the  spectral  density  of  a  given  (f,9)  particle  and  the  growth 
coefficient  ^  is  precomputed  in  WAVESET  for  5-knot  intervals  of  wind  speed.  Note 
also  that  ^  is  nonnegative  and  that  Sjn  is  evaluated  with  respect  to  the  local 
wind  direction  in  the  program. 

D.  NONLINEAR  TRANSFER 

The  nonlinear  transfer  is  represented  by  a  set  of  empirical  orthogonal 
functions  (eof's)  that  decompose  a  bofily  of  exact  calculations  of  the  transfer  for  a 
family  of  single-peaked  spectra  with  different  peakedness  and  angular  spread  s. 
The  theoretical  calculations  were  made  by  Hasselmann  and  Hasselmann  (1981).  The 
eof’s  were  constructd  by  T.  Barnett  at  Scripps  Institution  of  Oceanography.  The 


formula  used  to  evaluate  Sni  at  model  grid  points  is: 

(£,0)  *  D  [A(J/VR(«,SJ4  AST?  djr*  * 

An<  **,s  , 0bl'3] 

I 

where  f=  f/FM  and  ©i j  =  |9  -  0| .  The  following  parameters  ranges  are  used: 

0 *?  to-  ,  fisejj  *  iS# 

I  5  V  5  M  ,  6  If  *  o.i 

O  H>  $•  S  5  0,10  ,  4S  :  o  or 
o.(,L  5  £  ^  =  O  OJ 

The  mean  ABAR  and  standard  deviation  ASD  were  extracted  from  the  theoretical 
transfer  before  finding  eof's  so  that  all  parts  of  the  (J*,s)  range  would  have 
roughly  equal  weighting.  The  Bn  can  be  thought  of  as  the  eof's  and  the  An  as 
amplitudes  of  these  eof's. 

E.  LIMITING  SPECTRUM 

We  employ  the  well-known  Pierson-Moskowitz  spectrum  evaluated  at  the  local 
wind  frequency,  but  with  a  variable  saturation  range  parameter,  as  our  spectral 
limiter.  The  limiting  value  Ymax  for  given  f,9  is 

c>f  l'V+“4V]  <•.»*( 

TT  t 

The  logic  used  to  apply  the  limiter  is  described  in  Sec.  IV. C. 

IV.  SOLUTION  OF  THE  GROWTH  AND  PROPAGATION  EQUATIONS 
A.  OVERVIEW  OF  SUBROUTINE  INTRACT 

We  neglect  refraction  and  the  curvature  of  the  earth,  and  assume  that  the 
evolution  of  spectral  density  F  is  governed  by  a  radiative  transfer  equation  with 
the  general  form: 

ii.  F(*,y,+,e,t)  .  V'OF  '  i> 

dt 
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where  S  Is  the  net  source  function.  INTRACT  sets  up  the  source  terms  (mostly 
calculated  elsewhere)  that  apply  for  the  conditions  at  each  grid  point  (recall  that 
the  grid  Is  constant  for  all  directions  and  frequencies).  Then  the  ordinary 
differential  equation  given  above  is  Integrated  for  each  ray  point  (sampled  position 
along  one  of  a  family  of  parallel  characteristic  "curves"  for  the  given  direction 
and  frequency),  and  the  propagation  logic  is  applied. 

In  one  time  step  of  duration  At,  energy  of  a  given  frequency  can  move  a 
distance  Cg  &t,  where  Cg  is  the  group  velocity  of  waves  of  that  frequency.  In 
some  cases  Cg&t  may  be  small  compared  to  other  distances  of  interest.  The 
program  has  the  capability  of  holding  up  propagation  for  MXSTEP  (a  program  variable) 
steps,  while  still  allowing  energy  growth.  Then  consecutive  ray  points  are  (MXSTEP) 
Cg  &t  apart;  MXSTEP  may  be  different  for  each  frequency. 

The  energy  content  of  ray  points  on  the  incoming  boundary  is  fixed  at  a  very 
small  value  (<10-200).  When  MXSTEP>1,  this  condition  holds  until  the  energy 
packet  steps  to  the  next  ray  point.  Then  and  only  then  does  the  growth  equation 
apply.  Points  that  reach  the  outgoing  boundary  are  absorbed  perfectly;  there  is  no 
reflection  or  other  effect. 

B.  EQUATIONS  IMPLEMENTED 

Propagation  and  growth  are  logically  separated.  The  propagation  condition 
is 


X  -  Cgt  =  0 

where  x  stands  for  the  position  along  one  of  the  characteristic  rays.  As  described 
above,  propagation  takes  place  every  MXSTEP  time  steps  (MXSTEP  is  an  integer  1), 
so  the  propagation  condition  becomes 
&  x  =  Cg  (MXSTEP  .  &t) 


where  x  is  the  distance  moved.  As  Implemented  in  the  program,  the  array  that  holds 
energy  values  at  positions  along  rays  has  consecutive  elements  that  are  assumed  to 
be  AX  apart  along  a  given  ray;  the  first  point  on  each  ray  is  always  assumed  to  be 
on  the  incoming  boundary.  On  propagation  steps  an  (f,9)  particle  or  bundle  of 
energy  is  simply  moved  to  the  next-indexed  array  element.  (Implementation  detail: 
An  entire  family  of  rays  is  treated  as  one  array.  When  an  energy  packet  moves  from 
the  last  position  on  one  ray  to  the  next  array  location,  which  represents  the 
incoming  boundary  point  on  the  next  ray,  the  no-reflection  boundary  condition  is 
automatically  applied.) 

As  described  in  Section  III,  we  assume  that  the  net  source  function  S  is 
dominated  by  the  two  physical  processes: 

S  =  Sni  +  Sin 

where  $ni  is  synthesized  from  empirical  orthogonal  functions,  and 
Sm=((V4)Y 

where  Y  is  shorthand  for  the  directional  spectrum  at  the  particular  frequency, 
direction,  location,  and  time  under  consideration.  A  hybrid  integration  scheme  is 
used:  first-order  Euler  for  Sn]  and  effectively  an  implicit  trapezoidal 
integration  for  Sy^jnd-  The  composite  integraton  equation  for  advancing  to  time 
step  n+1  is 

1 * 

=  r  -  ♦  |  ,  -  ■']  r 

=  y"  *  s„,  dt  ♦  c  r*  y"*)  \  At 
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If  this  Is  a  propagation  step  then  Yn  and  Yn+1  are  evaluated  at  different 
locations.  The  quantities  (1  and  Sni  are  evaluated  at  the  same  location  as  Yn 
("growth  before  propagation").  The  limiting  form  of  the  spectrum  is  then  Invoked  as 
necessary  according  to  the  logic  discussed  in  the  following  section. 

C.  LOGIC  FOR  APPLYING  SPECTRAL  LIMITER 

The  maximum  spectral  density  at  any  ray  point  at  any  time  is  governed  by  the 
limiting  form  of  the  spectrum.  Numerous  logic  statements  are  used  to  decide  whether 
the  maximum  density  should  be  imposed.  These  statements  can  be  divided  for  purposes 
of  explanation  into  4  catagories:  subgrid  saturation  region,  swell  region, 
integrator,  and  saturation  region. 

1.  In  the  subgrid  region  all  other  calculations  are  by -passed  and  an 

( f ,0 )  particle  is  immediately  assigned  its  maximum  density.  Associated  with  each 
frequency  is  a  wind  speed  USAT  above  which  wave  growth  is  so  rapid  that  it  cannot  be 
resolved  in  the  model.  Measurements  of  net  growth  by  Barnett  and  Wilkenson  (1967) 
are  combined  with  the  notion  that  the  fractional  change  in  density  cannot  be  greater 
than  DSAT  in  one  time  step,  viz. 

&f/f  1  p(  M)„,t  **  * 

The  value  DSAT=1.5  was  used  in  the  earlier  version  of  the  model.  The  logic  is 
simply 

If  U>USAT( f ) ,  then  Y=Ymax 

2.  Swell  Region 

Wave  growth  below  the  wind  frequency  fu  is  prevented  if  the  spectrum 
is  already  greater  than  or  equal  to  the  limiter. 

If  f<fu(Xj )  and  Sn-|  ( )  >  0  and  Y>Ymax(xj ),  then  set 
^nl t xi )=0 ;  else  skip. 
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3.  Integrator 

The  incremental  change  DF  to  the  spectrum  during  a  time  step  is  forced 
to  be  nonnegative  (severely  notched  spectra  were  obtained  without  this  provision.) 
If  f>fu,  then  limit  DF  >0 
yn+l=  yn  +  DF 

4.  Saturation  Region 

Application  of  the  limiting  form  for  spectral  density  depends  on  fm, 
fu,  0  -  0u>  and  in  some  cases  the  density  at  the  previous  time  step. 

a.  If  f<fu(xi+1)  then  do  nothing 

b.  Else  (active  generation,  f >_fu ( x-j  +i ) 

(1)  If  IQ  -  9U|<90°  (forward  half-plane)  then 
If  Y<Ymax(xi+i)  then  do  nothing 

Else  (Y>Yrnax(xi+l) 

If  fffmlxi+i)  then 

If  Yn>ymax(x1>1)  then  set  Yn+l=yn 
Else  (Yn<Ymax(x1+1))  set  Y"+1=Y 

max(*  1+1) 

Else  (f>fm(x1+i))  set  Yn+l=Ymax(xi+i) 

(Zero  (infinite)  decay  time  for  components  above  (below)  local  peak) 

(2)  Else  | 6  -  Qu|>90°  (back  half -plane) 

If  f>fmUi+i)  then  set  Yn+I=  0  (saturation). 

Else  (f<fm(xi+i)}  do  nothing 

(3)  Eliminate  negative  values  -  Limit  Y  >_  lO'6^. 

D.  SOURCE  (Y)  AND  DESTINATION  (X)  CONCEPT  FOR  PROPAGATION 

The  program  uses  an  array  of  "spectral  words"  each  of  which  contains  the 
energy  density  (directional  spectrum)  at  a  certain  location  (ray  point),  for  the 
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frequency  and  direction  being  treated.  In  one  time  step  the  value  of  a  spectral 
word  generally  changes,  and  if  it  is  a  propagation  step  the  location  also  changes. 
The  program  uses  two  temporary  variables,  X  and  Y.  Roughly  speaking,  Y  is  a 
"source"  and  X  is  a  "destination".  As  each  point  in  the  array  is  processed  Y  starts 
off  with  the  value  for  location  xi  at  the  beginning  of  the  time  step.  Application 
of  the  net  source  function  to  the  old  spectrum  value  creates  X.  X  is  stored  into 
the  same  array  location  for  non-propagation  steps,  or  into  the  next  locations  for 
propagation  steps.  Using  previous  notation,  Y  =  Yn  and  X  =  Yn+1.  Y  is 
associated  with  location  x-j ,  X  is  associated  with  x-j+i  or  x^ ,  depending  on 
whether  propagation  does  or  does  not  take  place,  respectively. 

E.  TEMPORARY  CHANGES 

The  subgrid  saturation  logic  was  not  used  for  the  wave  model  intercompari son 
(Hasselmann  et  al . ,  1982).  We  made  two  temporary  changes  to  try  to  handle  problems 
that  arose  at  relatively  high  frequencies.  First,  the  saturation  condition  was 
applied  to  certain  (f,6)  particles  which  otherwise  would  not  receive  any  spectral 
density  (having  encountered  no  ray  point  'upwind'  where  DF  >  0).  The  logic  was 
if  propagation  step  and  f  fu(x-j)  and  f  >  fm(xi)  and  Y  <  10-20 
then  Y  =  Ymax 

Second,  an  exponential  filter  was  added  to  combat  a  sawtooth  effect  in  the  spectrum 
at  short  fetches.  The  recursive  filtering  scheme  was 
7n+l  =  Yn  +  DF  +  e_1/N  (X-Y) 

where  the  tildes  indicate  smoothed  values  and  N  =  MXSTEP/8  (where  MXSTEP  is  the 
number  of  time  steps  between  propagation  steps)  was  used. 

We  strongly  recommend  that  these  changes  be  removed  from  the  model,  and  for 
practical  calculations  replaced  with  a  subgrid  mechanism  such  as  described  in 
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IV. C.l.  These  temporary  changes  helped  very  little  and  actually  caused  minor 
problems. 

V.  USER  INSTRUCTIONS 

A.  WAVESET  data  cards 

[NOTE:  Variables  read  in  on  NAMELIST  data  cards  and  their  default  values 
(if  any)  are  defined  in  comments  at  the  beginning  of  WAVESET  (main  program).  The 
comments  show  the  NAMELIST  names  preceded  by  an  asterisk.  In  the  current  version 

the  asterisk  is  replaced  by  a  $  sign  in  column  2  (e.g.  _$SET,  $B0ND,  etc.).] 

The  first  card  has  *RUN  in  cols.  1-4.  Text  (title  information)  may  appear  in 
cols.  7-22. 

The  next  data  card  is  the  _$SET  card.  Important  parameters  are  listed  below: 
XOS,  YOS 
SZERO 
EZERO 
FQZERO 

SCALE  (=  .53996  for  distances  in  km) 

DTIME 

NSTEP  (=  0  normally) 

INITDT* 

Internal  grid  points  and  boundary  points  are  the  external  values  (from  data  cards) 
with  XOS  &  YOS  subtracted: 

<x.Y>int  =  (X,Y)ext  ‘  (XOS, YOS) 

♦Initial  date  and  time  may  be  specified  by  giving  ITIME  =  time  in  seconds  since 
midnight,  Jan.  1,  1900. 
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The  boundary  is  specified  by  the  _$B0ND  card.  For  a  single  basin  it  is 
sufficient  to  specify  NBOND  >  2  and  the  XBOND  &  YBOND  lists.  For  a  more  complicated 
situation  (e.g.  islands)  a  little  more  information  is  needed.  Example:  Suppose  the 
basin  is  10  x  10  with  lower  left  corner  at  (0,0),  containing  a  right  triangular 
island  extending  2  units  in  X  and  1  unit  in  Y,  with  the  right  angle  at  (2,4).  Then 
the  data  card  should  be 


$NB0ND  NBOND  =  7, 

XBOND  =  0,  0,  10,  10,  2,  2,  4, 

YBOND  =  0,  10,  10,  0,  4,  5,  4, 

NBSTART  (2)  =  5  SEND 

NBSTART  (1)  defaults  to  1  and  need  not  be  specified.  There  are  limits  of  200 
boundary  points  and  9  islands. 

Grid  points  are  specified  by  _$GRID  cards.  There  are  three  methods,  which  may 
be  combined: 

(1)  Specify  DX,  DY,  XS,  YS.  Then  points  will  be  generated  at  (XS  +  mDX,  YS  +  nDY ) , 
m  =  0,  +1,  +2,  — ;  n  =  0,  +1,  +2,  — .  All  points  that  fall  within  the  boundary 
will  be  retained.  The  program  counts  the  points  that  are  generated  and  retained. 

(2)  Give  LATTICE  4  0,  XLIST  and  YLIST  arrays,  and  NX  and  NY.  The  program 
generates  grid  points  at  all  intersections  (XLIST(n),  YLIST(m)),  n  =  1,2, — ,NX, 
m  =  1,2, — ,NY,  and  counts  the  points. 


20 


(3)  Explicitly  list  the  grid  points  (elements  of  XGRIO  and  YGRID  arrays).  To  use 
this  method  properly  it  is  necessary  to  know  NGRID  =  number  of  grid  points  that  have 
already  been  specified.  The  explicit  lists  begin  with  XGRID(NGRID  +  1)  and 
YGRID(NGRID  +1).  It  is  also  necessary  to  redefine  NGRID  to  account  for  the  number 
of  points  in  the  lists. 

Grid  point  input  is  terminated  when  the  last  _$GRID  card  leaves  LATTICE  and  at 
least  one  of  DX,DY  unset.  In  some  cases  an  extra  _$GRID  card  is  required  to  do 
this.  There  is  a  limit  of  500  grid  points. 

A  list  of  up  to  20  frequencies  may  be  read  in  using  a  _$FD  data  card.  The 
frequencies  must  be  ordered  monotonically  increasing.  The  remaining  cards  are  keyed 
to  the  frequency  list.  [A  mistake  in  the  original  program  that  prevented  the  20th 
frequency  from  being  read  was  corrected.] 

The  remaining  NAMELIST  cards  (_$JD  and  _$JK)  specify  the  families  of  rays.  The 
following  must  be  specified:  directions,  lateral  spacing  between  individual  rays, 
and  spacing  between  points  on  the  rays  (where  energy  is  localized  during 
propagation).  Families  of  parallel  rays  are  generated  in  directions  that  are 
derived  from  information  on  _$JD  card(s).  The  directions  are  usually  specified  for 
the  first  frequency  and  remain  specified  for  subsequent  frequencies,  although  they 
can  be  redefined.  The  first  nine  (9)  columns  of  the  card  have  a  mandatory  format: 

_ $ JD _ Q=j j ,  where  jj  is  the  index  of  the  specified  frequency  (e.g.,  0)  for  the  first 

frequency  in  the  _$FD  list).  Directions  are  specified  by  a  list  KDIR.  The  KDIR 
values  are  not  the  ray  directions.  Rather,  they  are  directions  to  cut  the  unit 
circle  into  "pie  slices".  The  ray  directions  are  along  lines  that  bisect  the  vertex 
angles  of  the  pie  slices.  This  guarantees  that  the  influence  of  a  given  direction 
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(in  integrals,  e.g.)  is  symmetric  about  that  direction.  Up  to  36  KDIR  values  may  be 
speci fied. 

The  value  MXSTEP  is  also  given  on  _$JD  cards;  it  determines  the  spacing  between 
points  on  a  ray.  A  natural  length  unit  for  the  propagation  problem  is  (group 
velocity)  x  (time  step).  Ray  points  may  be  placed  at  any  integral  multiple  of  this 
distance  (starting  on  the  incoming  boundary);  MXSTEP  is  the  multiplier.  The  same 
MXSTEP  value  remains  in  effect  until  it  is  redefined. 

So,  a  card  of  the  form 

_$JD_Q=0100,  KDIR  =  0,  30,  60,  90,  120,  330, 

MXSTEP  =  2  SEND 

specifies  the  following  situation:  beginning  with  the  first  frequency,  families  of 
parallel  rays  are  generated  in  the  directions  15°,  45°,  75°,  105°,  225°,  345°.  The 
spacing  between  successive  points  on  one  ray  is  twice  the  basic  natural  unit  of 
length.  These  parameters  remain  in  effect  until  they  are  redefined. 

As  diagrammed  in  1 1. A,  processing  effectively  involves  nested  loops:  for  each 
frequency  in  turn,  the  program  cycles  through  each  direction.  For  each  direction 
(for  each  frequency)  the  spacing  between  parallel  rays  may  be  specified,  using  _$JK 
cards.  The  first  eleven  (11)  columns  of  these  cards  must  have  the  format 
_$JK_Q=jjkk,  where  jj  gives  the  frequency  index,  as  for  _$JD  cards,  and  kk  gives  the 
direction  index.  There  are  three  ways  to  specify  ray  spacing. 

(1)  The  RAYSPAC  parameter  explicitly  gives  the  distance  between  rays.  In 
general,  the  number  of  rays  will  depend  on  their  direction. 

(2)  The  RAY  list  (up  to  100  values)  provides  a  method  of  obtaining  unequally 
spaced  rays,  among  other  things.  If  we  define  a  right-hand  coordinate  system  whose 
origin  coincides  with  the  origin  of  grid  and  boundary  coordinates,  and  whose 


positive  x-axis  extends  in  the  ray  direction,  then  the  RAY  list  gives  the 
y-coordinates  of  the  rays.  For  example,  in  the  45° 
direction  the  list 

RAY  =  -1,  0,  1,  3 
gives  the  rays  shown  at  right. 

(3)  The  NRD  parameter  defines  ray  spacing  as  a  fraction  of  the  transverse  width 
of  the  basin.  If  we  construct  a  pair  of  lines  in  the  ray  direction  that  are  tangent 
to  the  basin  at  the  maximum  possible  distance  on  either  side  of  the  origin;  then  the 
distance  between  these  lines  is  the  "transverse  width".  A  value  of  NRD  =  30  creates 
a  uniform  ray  spacing  1/30  of  this  value.  Then  the  distance  between  rays  will 
depend  on  direction,  but  the  number  of  rays  will  be  NRD  for  all  directions.  For 
example,  a  card  of  the  form 

_$JK_Q=0101,  RAYSPAC  =  25  SEND 

begins  with  the  first  frequency  and  the  first  direction,  and  specifies  rays  25 
external  length  units  (SCALE  x  nautical  miles)  apart.  This  specification  remains  in 
effect  for  all  frequencies  and  directions  until  it  is  redefined. 

It  is  permissible  to  intermix  _$JD  and  _$JK  cards  as  needed.  For  example, 

MXSTEP  may  be  redefined  at  higher  frequencies  to  prevent  the  ray  point  spacing  from 
becoming  significantly  less  than  the  grid  point  spacing  as  group  velocity  decreases. 

WAVESET  data  card  input  is  ended  by  a  card  of  the  form 
_$ST0P _ . 

B.  PROPAGR  data  cards 

The  first  data  card  has  _$RUN  in  columns  1-4,  with  text  in  columns  8-25. 

The  text  is  used  as  a  page  title  every  tenth  time  step. 
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The  next  card  is  a  _$DBTA  NAMELIST.  Only  two  of  the  variables  are  significant 
in  the  present  version: 

STEP  =  starting  time  step  (  =  0  normally  to  agree  with  WAVESET) 

STOP  =  stop  data  and  time  in  form  mm,dd,yy,hh,mm,ss 
Note:  One  step  is  performed  with  step  =  0,  i.e.  at  the  initial  time,  and  the 
stopping  condition  is  such  that  the  time  of  the  last  step  performed  must  be  greater 
than  the  stop  date-time.  So  the  number  of  processing  steps  is 

[((stop  date-time)  -  (start  date-time) )/(time  step)]  +  2 
The  next  data  card  has  the  form 
_$G0  text . 

The  preceeding  card,  and  the  complete  $DBTA  NAMELIST,  are  printed  before  processing 
begins. 

The  subroutine  that  produces  the  spectral  summary  tape,  and  associated  printout, 
is  controlled  by  a  _$SPRT  NAMELIST  card.  The  variables  are: 

NSP  =  Number  of  stations  (grid  points)  to  print  (_<  20) 

LOC  =  List  of  stations  to  print 
NPR  =  Time  step  interval  for  printing 

NSV  =  Time  step  interval  for  saving  data  on  spectral  suimary  tape 
For  the  purposes  of  this  output,  the  first  (dumniy)  time  step  is  treated  as  step  1. 
Also,  NPR  must  be  an  integral  multiple  of  NSV.  Thus,  if  DTIME  =  1200  (1200  seconds 
=  20  minutes;  value  set  in  WAVESET),  NSV  =  3,  and  NPR  =  24,  then  information  is 
written  on  the  spectral  summary  tape  every  hour  (after  3,  6,  9,  —  time  steps),  and 
printout  is  produced  every  8  hours. 

The  final  data  card  is 
_*ST0P _ , 

exactly  the  same  as  for  WAVESET. 
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C.  GRID/RAY  POINT  RELATIONSHIP 

For  each  frequency -direction  combination  there  are  effectively  two  grids  - 
the  “grid"  and  the  "rays".  The  former  is  used  to  (resolve  and)  specify  the  source 
term  fields,  and  to  provide  a  convenient  ouput  mechanism;  the  latter  forms  the 
spatial  mesh  on  which  the  discrete  approximation  to  the  radiative  transfer  equation 
is  solved,  so  it  must  resolve  the  wave  field.  The  grid  is  a  fixed  lattice  of 
points,  whereas  a  different  set  of  ray  points  is  used  for  each  frequency-direction 
combination. 

The  program  limits  the  number  of  grid  points  to  500.  Approximately  2000  ray 
points  in  each  set  are  allowed.  A  very  fine  ray  point  spacing  might  improve  the 
accuracy  of  the  integration,  except  that  the  source  terms  to  be  used  at  a  ray  point 
are  taken  to  be  those  calculated  at  the  nearest  grid  point.  Consequently,  there  is 
no  point  in  making  the  ray  point  spacing  much  smaller  than  the  grid  point  spacing. 
We  believe  that  a  good  rule  of  thumb  for  regular  grids  is  to  make  the  ray  point 
spacing  approximately  equal  to  the  nearest  neighbor  spacing  of  the  grid.  For 
example,  for  a  square  grid  this  is  about  2/3  of  the  grid  interval. 

Another  consideration  applies  an  upper  limit  to  the  spacing  of  points  along 
rays.  The  nonlinear  interaction  source  term,  calculated  at  grid  points,  involves 
integrals  over  the  spectral  density.  The  latter  is  calculated  at  ray  points.  So 
the  integrals  are  calculated  by  obtaining  the  directional  spectrum,  for  each 
direction  and  frequency  involved,  from  the  ray  point  nearest  to  the  grid  point  "at" 
which  the  calculation  is  performed.  This  mutual  relationship  --  energy  values  at 
ray  points  are  affected  by  source  terms  at  the  nearest  grid  point;  the  source  terms 
involve  spectral  values  at  nearby  ray  points  —  is  referred  to  as  the  "region  of 
influence"  concept. 
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Now  the  first  point  on  each  ray  is  on  the  incoming  boundary  and  is  constrained 
to  have  a  very  small  value.  Consequently,  it  is  important  to  prevent  the  grid 
points  closest  to  the  boundaries  from  "looking"  at  these  boundary  ray  points  in 
calculating  spectral  moments,  because  the  moments  would  be  underestimated  and  the 
effect  would  "feed  back"  to  other  ray  points  via  the  interaction. 

The  current  solution  to  this  problem  is  to  arrange  the  grid-ray  layout  so  that 
grid  points  close  to  boundaries  "look"  at  the  second  (or  later)  points  along  rays. 
WAVESET  could  be  modified  to  that  the  grid  point-ray  point  linkage  logic  would 
prohibit  a  grid  point  from  being  linked  to  a  boundary  ray  point.  The  modification 
is  not  trivial.  Inclusion  of  subgrid  scale  wave  growth  also  would  help  to  reduce 
this  problem. 

VI.  OUTPUT  FORMATS 

A.  CYCLE  AND  SUMMARY  TAPES 

In  each  time  step,  105-106  directional  spectrum  values  are  updated. 

The  source  term  calculation  and  interaction  logic  also  involves  thousands  of  other 
values.  This  extremely  large  number  of  variables  is  handled  by  the  use  of  external 
storage.  The  complete  state  of  the  system  is  kept  in  a  file,  which  is  designated 
the  cycle  tape.  Spectral  moments  are  also  kept  on  a  separate  summary  tape.  There 
are  two  cycle  tapes,  A  and  B,  and  two  summary  tapes,  C  and  D.  In  one  cycle  step 
PROPAGR  advances  the  solution  one  time  step  by  processing  A  and  C  to  produce  B  and 
D,  containing  the  updated  values.  In  the  next  time  step  the  roles  are  reversed. 

The  initial  state  is  set  up  by  WAVESET  (cycle  tape  only;  the  initial  sunroary  tape  is 
created  by  PROPAGR's  initialization  routine).  The  cycle  and  summary  tapes  at  the 
end  of  processing  serve  as  the  primary  means  of  output  for  the  final  state. 


26 


The  cycle  tape  format  is  shown  in  Table  3.  There  is  one  group  of  JREC,  JKREC, 
JSREC  records  for  each  frequency;  all  of  the  groups  precede  the  SREC  record.  The 
JREC,  JSREC,  and  SREC  records  also  appear  on  the  summary  tape.  A  detailed 
description  of  ecah  type  of  record  on  the  cycle  tape  is  given  in  Table  4. 

B.  SPECTRAL  SUMMARY  TAPE 

The  spectral  summary  tape  contains  the  directional  spectrum  and  a  few  other 
pieces  of  information.  Unlike  the  cycle  and  summary  tapes,  the  spectral  summary  is 
cumulative  --  information  is  written  to  it  every  NSV  time  steps  (V.B).  At  each  time 
step  for  which  information  is  written  to  tape,  a  5-word  ID  record  and  several 
spectral  records  —  one  for  each  direction  for  each  frequency  —  are  written.  For 
convenience  in  interpretation,  spectra  are  evaluated  at  (all)  grid  points,  using 
nearest-neighbor  interpolation.  The  structure  of  the  individual  records  is 
described  in  Table  5. 

VII.  SUMMARY 

The  present  model,  a  modified  version  of  the  original  model  by  Barnett  et  al . 
(1969),  received  its  first  workout  in  the  wave  model  intercomparison  (Hasselmann  et 
al.,  1982).  The  model  imitated  the  JONSWAP  laws,  although  it  is  not  based  on  them 
directly.  The  model  produced  good  results  for  many  ideal  cases  of  wave  generation 
from  a  simple  fetch-limited  situation  to  the  complex  situation  of  a  rapidly  moving 
hurricane.  We  found  a  lot  of  chatter  in  the  spectral  values  (and  the  various 
moments)  at  fetches  less  than  about  two  grid  lengths  because  the  gradients  in  the 
wave  field  were  not  resolved.  A  short  test  with  a  5-fold  increase  in  spatial 
resolution  cured  the  problem  but  was,  of  course,  impractical.  We  also  caused  some 
isolated,  undesirable  peaks  to  appear  in  complex  situations  through  the  temporary 
changes  discussed  in  IV.  More  than  likely  both  of  these  problems  would  be  cured  by 
using  some  mechanism  for  subgrid  wave  growth  as  suggested  in  IV.C.l.  (All  in  all 
the  view  we  took  for  the  intercomparison  was  too  pristine.) 


TABLE  3 


CYCLE  TAPE  FORMAT 


RECORD 

REMARKS 

LENGTH 

GRIDREC 

Mi  sc.  data  &  grid  point  coordinates 

2  x  NGRID  +  9 

BONDREC 

Boundary  description 

2  X  NBOND  +  13 

GREC 

Geometry  description  used  to  set  up 
linkages 

LGREC  <  2000 

JREC 

J  (frequency)  coefficients 

LJREC  (=  812) 

JKREC 

Spectral  densities,  one  record  for 
each  direction 

<  2049 

JSREC 

Sums  over  directions  for  this  frequency 

NGRID  +  2 

SREC 

Sums  over  frequency  and  direction 

3  x  NGRID  +  12 
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TABLE  4 


CYCLE  TAPE  RECORDS* 


GRIDREC 

0 

1 

2 

3 

4 

5 

6 

7 

8 

NG+8 

2NG+8 


ID=200g 

NGRID 

SCALE 

CDIST 

DTIME 

ITIME 

XOS 

YOS 

XGRID  (NGRID)  i 
YGRID  (NGRID)  ) 
ID 


BONDREC 

0  ID=201g 

1  NBOND 

2  NSTART  (10) 

12  XBOND  (NBOND) 

NBOND+12  YBOND  (NBOND) 

2NB0ND+12  ID 


(sec) 

(sec) 


INTERNAL  COORDINATES 


INTERNAL  COORDINATES 


GREC 

0  ID=400q  for  1st  grid  point  |  NN  I  MORE+  I  BYTE  I  BYTE  I 

'  15  '  TT  15  T5 - 

1  GREC  (LGREC) 

for  NGRID  1  NN  I  MORE*  |  BYTE  1  BYTE  1 

♦Note  that  word  counts  in  records  start  with  0,  not  1. 

+Pointer  to  continuation  word. 
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TABLE  4  (Cont.) 


GREC  (MORE)  |  BYTE  BYTE  1  BYTE  |  BYTE 

LGREC+1  ID 

'  • 

» 

JREC 

0  ID=600+J 

ft 

I  WINDMAXj  IWMAX,  RSTA,  EIGHT(8)  wind  above  which 

saturation  occurs,  saturation  value, 
extra  stuff 

II  ALPHA  (40,  10)-« — interaction  coefficients 

411  BETA  (40,  10)  -• — for  this  scheme. 

811  ID 

JKREC 

0  ID  =  |  0  |  NSTEP  |  J  |  K  1  One  per  freq.,  per  direction 

15  15  15  15 


i.e..  Lots! 

1 

NP 

No.  of  points 

2 

NNB 

Number  of  Neighbor  Words 

3 

MS 

Step  counter  for  this  J,  K 

modulo  MXSTEP 

4 

MXSTEP 

5 

KDIR 

0  to  360  integer 

6 

DDIR 

AKDIR/360  (sum  to  1.00) 

7 

DFREQ 

8 

FREQ 

9 

10 

X  I 

>  unused 

X  ) 
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TABLE  4  (Cont.) 


11  SPECT(NP) 


NP  +  11  NBT(NNB) 


NP+NNB+11  ID 

OSREC 

0  ID  = 

1  SEJ(NGRID) 

NGRID+1  ID 


L 

B 

I 

0 

1 

N 

1  |  N  i 

1 

ENERGY 

K 

0  I  D  I 

LX=1, NGRID  I 

42 

3 

15 

18 


LX  Index  of  nearest  grid  point 

LINK=1  this  ray  pt  also  nearest  to  that 

grid  point,  else  LINK=0 

B0ND=1  this  is  a  boundary  point  (incoming) 

Neighbor  table.  For  each  grid  point 

LX  that  is  not  the  nearest  neighbor  of 

any  ray  point,  i.e.  a  grid  pt  w/no  links, 

this  gives  the  ray  point  IX  that  is 

nearest. 

NBT  I  IX  |  LX  | 

30  30 

IX  closest  to  LX  where  LX  not  closest  to  IX. 

|  0  |  NSTEP  |  J  |  100  | 

(Integrals  of  energy  over  direction  but 
'for  this  freq. 


SREC 

0 

1 

2 

3 

4 

5 


JO  =  NSTEP*230+100B 

NGRID 

JMAX 

EZERO 

FQZERO 
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TABLE  4  (Cont.) 


6  THZERO 

7  I TIME 

8  X  \ 

9  X  )  not  used 

10  X  ) 

11  SE  (NGRID) 

NG+11  SFQ  (NGRID) 

2*NG+11  STH  (NGRID) 

3*NG+11  ID 

[NG  =  NGRID] 
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TABLE  5 


SPECTRAL  SUMMARY  TAPE  DESCRIPTION 

ID  RECORD 


1 

ID 

7777777777777777STEP  (octal) 

2 

LEN 

total  length  (=5) 

3 

ITIME 

time  in  seconds  since  1/1/00 

4 

JMAX 

number  of  frequencies 

5 

NGRID 

number  of  grid  points 

DIRECTIONAL  SPECTRUM  RECORD 

1 

ID 

10008  x  JFREQ  +  KDIREC 

2 

LEN 

total  length  =  4  +  [(NGRID  +  l)/2] 

3 

DF,  F 

Af,  f  (30  bits  each,  packed) 

4 

DTH,  K 

2*a6,  9*  (30  bits  each,  packed) 

5 

SP1,  SP2 

Spectrum  at  first  two  grid  points 

(30  bits  each,  packed) 

LEN  •  .  •  SPN  spectrum  at  last  (one  or  two)  grid  point(s) 


♦Angle  in  degrees. 
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An  open  question  is:  how  should  the  model  respond  to  a  rapidly  changing  wind 
direction?  The  present  version  probably  produces  too  large  a  drop  in  the  variance 
of  the  spectrum  when  there  is  a  change  in  the  wind  field,  although  the  rate  of 
turning  of  the  spectrum  (directional  relaxation)  hit  the  middle  ground  in  the 
intercomparison  study.  The  logic  described  in  Sec.  IV. c  needs  to  be  modified. 
Instead  of  using  the  local  peak  frequency  (and  the  wind  frequency)  to  define  a 
frequency  border  for  applying  the  spectral  limiter,  we  believe  it  would  be  worth 
trying  a  border  that  depends  on  the  average  wave  direction,  namely,  ^=0.13 
g/(U cos(9u~0)  for  |9u-9|<90°.  Comparison  of  model  behavior  with 
observations  (and  limited  theory)  is  necessary  here. 

The  real  test  of  the  model,  we  believe,  would  come  in  an  actual  hindcasting 
study.  Application  of  rigorous  techniques  for  comparing  model  results  and  data 
would  be  very  helpful.  A  key  part  of  such  a  study  would  be  trying  to  assign  how 
much  of  the  ‘error1  might  be  due  to  model  inadequacies  vs.  how  much  might  be  due  to 
uncertainties  in  the  wind  field. 
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